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Abstract. We present an extension of our GPGCD method, an iterative 
method for calculating approximate greatest common divisor (GCD) of 
univariate polynomials, to multiple polynomial inputs. For a given pair 
■^^ ■ of polynomials and a degree, our algorithm finds a pair of polynomials 

,— n ' which has a GCD of the given degree and whose coefficients are perturbed 

from those in the original inputs, making the perturbations as small as 
possible, along with the GCD. In our GPGCD method, the problem of 
approximate GCD is transferred to a constrained minimization problem, 
then solved with the so-called modified Newton method, which is a gen- 
eralization of the gradient-projection method, by searching the solution 
iteratively. In this paper, we extend our method to accept more than two 
polynomials with the real coefficients as an input. 
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1 Introduction 



f^ ' For algebraic computations on polynomials and matrices, approximate algebraic 

algorithms are attracting broad range of attentions recently. These algorithms 
take inputs with some "noise" such as polynomials with floating-point number 
coefficients with rounding errors, or more practical errors such as measurement 
errors, then, with minimal changes on the inputs, seek a meaningful answer 

'V^j . that reflect desired property of the input, such as a common factor of a given 

degree. By this characteristic, approximate algebraic algorithms are expected to 
be applicable to more wide range of problems, especially those to which exact 
algebraic algorithms were not applicable. 

As an approximate algebraic algorithm, we consider calculating the approxi- 
mate greatest common divisor (GCD) of univariate polynomials, such that, for a 
given pair of polynomials and a degree d, finding a pair of polynomials which has 
a GCD of degree d and whose coefficients are perturbations from those in the 
original inputs, with making the perturbations as small as possible, along with 
the GCD. This problem has been extensively studied with various approaches 
including the Euclidean method on the polynomial remainder sequence (PRS) 
(PQ, [2]> [3]) j the singular value decomposition (SVD) of the Sylvester matrix (g], 



[5J), the QR factorization of the Sylvester matrix or its displacements ([5J, [7], 
[8 ]), Pade approximation [5], optimization strategies ([IS], [II] , [H], [13], [II])- 
Furthermore, stable methods for ill-conditioned problems have been discussed 

m, na, csd- 

Among methods in the above, we focus our attention on optimization strate- 
gies. Already proposed algorithms utilize iterative methods including the Levenberg 
Marquardt method [10] , the Gauss- Newton method [14] and the structured total 
least norm (STLN) method ([11], [12] )■ Among them, STLN-based methods have 
shown good performance calculating approximate GCD with sufficiently small 
perturbations efficiently. 

In this paper, we discuss an extension of the GPGCD method, proposed by 
the present author (|17j. |21)~). an iterative method with transferring the original 
approximate GCD problem into a constrained optimization problem, then solv- 
ing it by the so-called modified Newton method [18] , which is a generalization 
of the gradient-projection method [19] . In the previous papers ([17], [Hj), we 
have shown that our method calculates approximate GCD with perturbations 
as small as those calculated by the STLN-based methods and with significantly 
better efficiency than theirs. While our previous methods accept two polynomials 
with the real or the complex coefficients as inputs and outputs, respectively, we 
extend it to handle more than two polynomial inputs with the real coefficients 
in this paper. 

The rest part of the paper is organized as follows. In Section [2] wc transform 
the approximate GCD problem into a constrained minimization problem for the 
case with the complex coefficients. In Section [3J we show details for calculating 
the approximate GCD, with discussing issues in minimizations. In Section^ we 
demonstrate performance of our algorithm with experiments. 



2 Formulation of the Approximate GCD Problem 

Let P\{x), . . . , P n (x) be real univariate polynomials of degree d±, . . . , d n , respec- 
tively, given as 

P A X ) =Pd! x H Pi x + Po > 

for i — 1, . . . , n, with minjdi, . . . , d n } > 0. We permit Pi and Pj be relatively 
prime for any i ^ j in general. For a given integer d satisfying minjdi, . . . , d n } > 
d > 0, let us calculate a deformation of Pi (x), ... , P n (x) in the form of 

Pi(x) = Pi(x) + APi(x) = H{x) ■ Pi(x), 

where APi(x) is a real polynomial whose degrees do not exceed di, respec- 
tively, H(x) is a polynomial of degree d, and Pi(x) and Pj(x) are pairwise rel- 
atively prime for any i ^ j. In this situation, H (x) is an approximate GCD of 
P\{x), . . . ,P n (x). For a given d, we try to minimize ||APi(x)||2 + - • ■+||Z\P rl (a;)||2, 
the norm of the deformations. 



For a real univariate polynomial P(x) represented as P(x) = p n x n +- ■ -+pqx°, 
let Ck(P) be a real (n + k, k + 1) matrix defined as 



C k (P) = 



(Vn \ 



pa Pn 



V PoJ 

fc+1 
and let p be the coefficient vector of P(x) denned as 

P= {Pn,---,Po)- 



(1) 



In this paper, for a generalized Sylvester matrix, we use a formulation by 
Rupprecht [20l Sect. 3]. Then, a generalized Sylvester matrix for Pi,...,P n 
becomes as 



N(P x ,...,P n ) = 



(C dl ^{P 2 ) C d2 _i(Pi) 
C dl _!(P 3 ) C^^PO 



VC dl _i(P n ) 







\ 





C d „_i(Pi)/ 



, (2) 



and the fc-th subresultant matrix (with min{(ii, . . . , d n } > fc > 0) is also defined 
similarly as 



N k (P 1 ,...,P n ) 

/C dl _i_ fe (P 2 ) Ck_i_ k (Pi) o 

C dl _i_ fe (P 3 ) Cd,_i_ fc (Pi) 








\ 



with 
rows and 



\C dl _i_ fe (P n ) ••• c d „_ 1 _ fc (Pi)y 

n- = ^i + c?2 + • • • + d n - (n - I)fc + (n— 2)di 
Ck = di + d 2 + ■ ■ ■ + d„ — n-k 



(3) 

(4) 
(5) 



columns. 

Calculation of GCD is based on the following fact. 

Proposition 1 (Rupprecht [20, Proposition 3.1]). JV/.(Pi, . . . ,P„) has full 
rank if and only if deg(gcd(Pi, . . . , P n )) < k. 

Thus, for a given degree d, if Nd-i(Pi, ■ ■ ■ ,P n ) is rank-deficient, then there 
exist real univariate polynomials Ui(x), . . . , U n (x) of degree at most d\—d, . . . , d n - 
d, respectively, satisfying 

UxPi + U Z P! = 0, (6) 



for i = 2, . . . , n. In such a case, if Ui and Uj are pairwise relatively prime for 
any i ^ j, then P = -^ = — S. = • • . = — £1 becomes the expected GCD. 
Therefore, for given polynomials Pi, . . . , P„ and a degree (i, our problem is to 
find perturbations Z\Pi , . . . , AP n along with cofactors U\ , . . . , U n satisfying ^ 
with making ||^Pi(x)||| + • • • + ||AP n (x)||| as small as possible. 
By representing Pj(x) and Ui{x) as 



U t (x) = u£_ d a;*-" + • • • + ufx + uf 



(7) 

we express the objective function and the constraint as follows. For the objective 
function, ||APi(x)||2 + • • • + ||AP„(a;)||! becomes as 



n [ di 

\\AP 1 (x)\\l + ■■■ + \\AP n (x)g = E E (p? -P?) ) ■ ( 8 ) 

i=l [j=0 

For the constraint, ((6]) becomes as 

iV d _i(P 1 ,---,Pn)- t (wi,---,Wn) = 0, (9) 

where Ui is the coefficient vector of Ui(x) defined as in (TIJ. Furthermore, we add 
another constraint for the coefficient of Ui (x) as 

ra! + - + HEy§ = i, (io) 

which can be represented together with (j9|) as 

GwW.fo 1 )-''* "-"^ <"» 

where (fT0|) has been put on the top of (|9|). Note that, in (fTT|) . we have total of 

d = di + --' + dn-(n-l)(d-l) + (n- 2)di + 1 (12) 

equations in the coefficients of polynomials in ([7|) as a constraint, with the j-th 
row of which is expressed as gj = 0. 
Now, we substitute the variables 

'd 1 ) d 1 ) d n ) d n ) .,« ,.« „» „,(") 



KPdi ' m ■ m ' fo 1 ■ ■ ■ 'Pd„ ' • • ->^0 ' u di-d' ' • ■ ' u ' • ■ • ' u d„-d' • ■ • ' u /' \ LO ) 

as x = (xx, . . . ,£2(dH hd„)+(2-d)n)) then ([8|) and (fTT|) become as 

/( 3; ) = (x 1 -p( 1 i ) ) 2 + ... + (.T dl -rf ) ) 2 + --- 

h (x dl+ ... +dn _ 1+n - p d n ^>) 2 H h (a;d 1 +...+d„_ 1 +d„+n - Po n) ) 2 > (14) 

g(x) = t (g 1 {x),...,g i (x)) = 0, (15) 

respectively, where J in (TT5|) is defined as in (TT2"1) . Therefore, the problem of 
finding an approximate GCD can be formulated as a constrained minimization 
problem of finding a minimizer of the objective function f(x) in (Tl4l) . subject to 
g(x) =0 in flTSJ}. 



3 The Algorithm for Approximate GCD 

We calculate an approximate GCD by solving the constrained minimization 
problem (fH|) . (IT5|) with the gradient projection method by Rosen [TH] (whose ini- 
tials become the name of our GPGCD method) or the modified Newton method 
by Tanabe [TB] (for review, see the author's previous paper [T7]). Our preceding 
experiments ([T7] Sect. 5.1], [5TJ Sect. 4]) have shown that the modified Newton 
method was more efficient than the original gradient projection method while 
the both methods have shown almost the same convergence property, thus we 
adopt the modified Newton method in this paper. 

In applying the modified Newton method to the approximate GCD problem, 
we discuss issues in the construction of the algorithm in detail, such as 

— Representation of the Jacobian matrix J g {x) and certifying that J g (x) has 
full rank (Sect. EH]), 

— Setting the initial values (Sect. [32]), 

— Regarding the minimization problem as the minimum distance problem 
(Sect. E3J|, 

— Calculating the actual GCD and correcting the coefficients of P, (Sect. [ 

as follows. 



3.1 Representation and the rank of the Jacobian Matrix 

By the definition of the constraint (1151) . we have the Jacobian matrix J g {x) (with 
the original notation of variables (Tl3]) for a;) as 

/ ■•• 

C dl (U 2 ) C d2 (U t ) ■•• 
J a (x)= C *dUs) C d3 (U x ) 

\C dl (U n ) ••• C dn (U x ) 

2- t u x 2-*m 2 2-'m 3 ••• 2- t u n \ 
C dl ^ d {P 2 ) C d2 _ d {P x ) • • • 

C dl _ d (P 3 ) C d3 - d (P x ) 

C dl - d (P n ) ••• C dn _ d (P x )J 

which can easily be constructed in every iteration. Note that the number of rows 
in J g {x) is equal to d in (fT2"|) . which is equal to the number of constraints, while 
the number of columns is equal to 2(d x + • • • + d„) + (2 — d)n, which is equal to 
the number of variables (see ([15])). 

In executing iterations, we need to keep that J g {x) has full rank: otherwise, 
we are unable to decide proper search direction. For this requirement, we have 
the following observations. 



Proposition 2. Assume that we have degd < min{di, . . . , d n } — 1 and dcgUi > 
1 for i = l,...,n. Let x* G V g be any feasible point satisfying (|15[) . Then, if 
the corresponding polynomials do not have a GCD whose degree exceeds d, then 
J g {x*) has full rank. 



Proof. Let 



■d 1 ) 



d 1 ) 



An) 



d n ) J 1 ) 



(1) 



,W 



.(n). 



* f -v 1 ) -k 1 ) ~v l J ~\ !i ) \ L ) y 1 ) \ rl ) \ tl )\ 

x — \Pdx '• --iPo '■■■iPd„i---'P0 > w di-d> • ■ • ' M ! ■ • ■ : u d„-d' ■ • ■ ' u J 

as in (1131) , with its polynomial representation expressed as in (JT]) (note that this 
assumption permits the polynomials Pi{x) to be relatively prime in general). To 



verify our claim, we show that we have rank( J g {x*)) = d = d\ + • 



(n- 



l)(d— l) + (n — 2)di + 1 (see (fT2|) V Let us divide J g (x*) into two column blocks 
such that J g {x*~) = (Jl | Jr.), where Jl and Jr are expressed as 



Jl 



( ■•• 

C dl ([/ 2 ) C d2 (C/i) ••• 
c dl (u 3 ) C d3 (C/i) 



\ 



Jr = 



\C dl (U n ) ••• C dn {Ux)J 

( 2 • 'mi 2 • *m 2 2 • *m 3 • • • 2 • *«„ \ 

c dl _ d (P 2 ) c^-dtPi) o ••• o 

C dl _ d (P 3 ) C d3 _ d (Pi) 



\Cdx-«j(^n) ••• C d „- d (Pl)J 

respectively. Then, we have the following lemma. 

Lemma 1. We have rank(JL) = d = di H \-d n — (n — l)(d — 1) + (n — 2)di. 

Proof. Let J be a submatrix of Jl by eliminating the top row. Since the number 

of rows in Jl is equal to d = di H h d n — (n — l)(d — 1) + (n — 2)di, we show 

that J has full rank. 

For i = 2, . . . , n, let us divide column blocks C^ (C/i) and C^ (C/i) as 

d+1 di-d 

C dl (C/ 4 )= (C dl (Ui) h C dl (Ui) K ) , 

C A (J7i)= (^([/Ol^^Or) 
'C d (tfi)\ 



C dl (C/i)i 



\ }d+i 

C dl _ d -l(C/i) / }di+*-2rf 



(16) 



Cd^Uih 







}di-d 



C d ,(t/: 



ljR 



' \ }d+i 

.Gdk-d-xiUi)) }d!+d z -2d 



(17) 



respectively, thus J is expressed as 



J = 



f C dl (U 2 ) h C dl (U 2 )n 
C dl (U 3 )L C dl (U 3 ) R 

\C dl (U n )i, C dl (U n ) R 



C<2 2 (/7i)l Cd 2 {Ui)s. 







C < i a (l7i)LC d8 (17i)L 



• 









\ 





CdAUihCdjU^nJ 



Then, by exchanges of columns, we can transform J to J = (Jl Jr) , where 



Jl = 



Jr = 



fC dl (U 2 )-L C d2 ([/!) L ••■ \ 
C dl (t/ 3 ) L C d3 (C/i) L 



\C , dl (C/„) L ■■• C dn (Ui) h J 

/Cd 1 (U2)nC d2 (U 1 ) R •■■ \ 

C dl (C/ 3 )R C d3 (C7i) R o 



\C dl (U n )j 







CdjU^nJ 



We see that nonempty rows in J R consist of N(Ui, ■ ■ ■ , C/ n ), a generalized 
Sylvester matrix for Ui,...,U n (see ©). By the assumption, Ui,...,U n are 
pairwise relatively prime, thus, by Prop. [TJ rank( Jr) is equal to the number of 
nonempty rows in Jr, which is equal to d 2 + ■ ■ ■ + d n + (n — l)(di — 2d) (see (J16I) 
and CEZ])). 

On the other hand, in Jl, column blocks Cd 2 {Ui)h,Cd 3 (Ui)i,, . . . , Cd n {Ui)i, 
are lower triangular matrices with d + 1 diagonal elements, which shows that 
rank(Ji) is equal to the sum of the number of columns in C<2 2 ([/i)l, C^CAJl, ■ • ■, 
Cd n (Ui)l, which is equal to (n — l)(d + 1). 

Furthermore, we see that the row position of diagonal elements in Cd 2 {Ui)i,, 
Cd 3 (Ui)-L, . . . , Cd n (Ui)i, correspond to the position of the empty rows in Jr, 
thus the columns in Cd 2 (V r i)j J , Cd 3 {U{)i^ . . . , C c ;„(t r i)L are linearly independent 
along with the columns in Jr. Therefore, we have 

rank(J) = rank( Jl) + rank(Jj{) = d\ -\ 1- d n — (n — l)(d — 1) + (n — 2)di, 

which proves the lemma. 

Proof of Proposition [Z\ (continued). By the assumptions, we have at 
least one nonzero coordinate in the top row in Jr, while we have no nonzero 
coordinate in the top row in Jl, thus we have rank(J g (a;)) = d\ + ■ • • + d n — 
(n — l)(d — 1) + (n — 2)d\ + 1, which proves the proposition. □ 



Proposition^ says that, under certain conditions, so long as the search direc- 
tion in the minimization problem satisfies that corresponding polynomials have 
a GCD of degree not exceeding d, then J g {x) has full rank, thus we can safely 
calculate the next search direction for approximate GCD. 

3.2 Setting the Initial Values 

At the beginning of iterations, we give the initial value Xq by using the singu- 
lar value decomposition (SVD) [22] of N d -i(Pi, . . . ,P n ) (see ©) as N d -i = 
U S f V,U = (wi,...,w Cd _ 1 ), S = diag(ci,.. . jO-^J, V = (v 1 , . . . , v Cd _ ± ), 
where Wj G R r<i_I , Vj G R Cd_1 with r^ and Cfc as in (JH) and ([5]), respectively, 
and S = diag(<7i, . . . , cr Cdl ) denotes the diagonal matrix with the j-th di- 
agonal element of which is <jj. Note that U and V are orthogonal matrices. 
Then, by a property of the SVD [22] Theorem 3.3], the smallest singular value 
<r Cdl gives the minimum distance of the image of the unit sphere S c<i_1_ , given 
as S^- 1 " 1 ={i G R^- 1 | ||jc|| 2 = 1}, by N d -i(P x ,...,P n ), represented as 
N d -i ■ S "- 1 ^ 1 = {N d -ix | x e R^- 1 , ||jc|| 2 = 1}, from the origin, along with 
(7 Cdl w Cdl as its coordinates. Thus, we have N d -i ■ v Cdl = a Cdl w Cdl . For 

V, = ^SU-.^.-.^U-.^), let Ui(x) = u®_ da *- d + -..+ 
Uq x° fori = 1, . . .,n. Then, U\{x),. . ., U n (x) give the least norm of UiPi + UiPi 

satisfying \\Ui\\\ H h ||f/ n ||| = 1 by putting Ui{x) = Ui{x) in Q. 

Therefore, we admit the coefficients oi Pi, . . . ,P n , Ui, ■ ■ ■ , U n as the initial 
values of the iterations as 

_ I (1) (1) (n) (n) -(1) -(1) _(n) _(n)x 

X — \P dl J • • • J Pq j • • • iPd„ ) • • • iPo i u di-d' ■ • ■ ' u ) ■ • ■ i u d n -di ■ • ■ ' M )• 

3.3 Regarding the Minimization Problem as the Minimum Distance 
(Least Squares) Problem 

Since we have the object function / as in (fl4|) . we have 

Vf(x) = 2- t (x 1 - Pd \\...,Xd 1 -p ( o\--., 

(") (™) n rA 

£di + ---+d„_i+n - P d?i , • ■ • , aJd 1 +..-+d„_ 1 +d„+n — Pq , U, . . • , UJ. 

However, we can regard our problem as finding a point cc G V^ which has the min- 
imum distance to the initial point Xq with respect to the (xi, ... , x dl -\ \- dn _ x + dn +n)- 

coordinates which correspond to the coefficients in Pi(x). Therefore, as in the 
case for two polynomials (see the author's previous papers ([12], [21])), we change 
the objective function as f(x) = ^f(x), then solve the minimization problem of 
f(x), subject to g(x) = 0. 

3.4 Calculating the Actual GCD and Correcting the Deformed 
Polynomials 

After successful end of the iterations, we obtain the coefficients of Pi{x) and 
Ui{x) satisfying (j6]) with Ui(x) are relatively prime. Then, we need to compute 



the actual GCD H(x) of Pi(x). Although H can be calculated as the quotient of 
Pi divided by Ui, naive polynomial division may cause numerical errors in the 
coefficient. Thus, we calculate the coefficients of H by the so-called least squares 
division [13] , followed by correcting the coefficients in Pi by using the calculated 
H, as follows. 

For polynomials Pi, and Ui represented as in and H represented as 

H(x) = h d x d H h h x°, 

solve the equations HUi = Pi with respect to H as solving the least squares 
problems of a linear system 

C d (U i ) t (h d ...,h ) = t (p%,---,p { o ) ). (18) 

Let Hi(x) £ R[x] be a candidate for the GCD whose coefficients are calculated 
as the least squares solutions of (|18l) . Then, for i = 2, . . . , n, calculate the norms 
of the residues as 



|2 



r i = 53||P i -Jf i 17 J -| 

.7=1 

and set the GCD H(x) be Hi(x) giving the minimum value of r^ so that the 
perturbed polynomials make the minimum amount of perturbations in total. 

Finally, for the chosen H(x), correct the coefficients of Pi(x) as Pi(x) = 
H(x) ■ Ui(x) for i = l,...,n. 

4 Experiments 

We have implemented our GPGCD method on the computer algebra system 
Maple and compared its performance with a method based on the structured 
total least norm (STLN) method [TT] for randomly generated polynomials with 
approximate GCD. The tests have been carried out on Intel Core2 Duo Mobile 
Processor T7400 (in Apple MacBook "Mid-2007" model) at 2.16 GHz with RAM 
2GB, under MacOS X 10.5. 

In the tests, we have generated random polynomials with GCD then added 
noise, as follows. First, we have generated a monic polynomial Pq(x) of degree 
to with the GCD of degree d. The GCD and the prime parts of degree m — d 
are generated as monic polynomials and with random coefficients c £ [—10, 10] 
of floating-point numbers. For noise, we have generated a polynomial Pn(x) of 
degree to — 1 with random coefficients as the same as for Pq(x). Then, we have 
defined a test polynomial P(x) as P(x) = Pq(x) + n P e , F Mi Pn(x), scaling the 
noise such that the 2-norm of the noise for P is equal to ep. In the present test, 
we set ep = 0.1. 

In this test, we have compared our implementation against a method based on 
the structured total least norm (STLN) method |11] , using their implementation 
(see Acknowledgments). In their STLN-based method, we have used the proce- 
dure R_con_mulpoly which calculates the approximate GCD of several polyno- 
mials in R[x]. The tests have been carried out on Maple 13 with Digits=15 



Table 1. Test results for (m,d,n): n input polynomials of degree m with the 
degree of approximate GCD d. See Section 2] for details. 



Ex. 


(m, d, n) 


#Fail 


Error 


^Iterations 


Time (sec.) 


STLN 


GPGCD 


STLN 


GPGCD 


STLN 


GPGCD 


STLN 


GPGCD 


1 


(10,5,3) 








2.31e-3 


2.38e-3 


5.50 


11.2 


1.17 


0.45 


2 


(10,5,5) 








5.27e-3 


5.22e-3 


4.70 


13.5 


3.10 


1.53 


3 


(10,5,10) 








5.48e-3 


5.62e-3 


4.40 


17.9 


12.49 


8.59 


4 


(20,10,3) 








5.17e-3 


5.40e-3 


4.50 


12.0 


3.35 


1.52 


5 


(20,10,5) 








5.89e-3 


5.85e-3 


4.40 


12.7 


10.37 


4.97 


6 


(20, 10, 10) 





1 


6.31e-3 


6.20e-3 


4.00 


25.6 


44.62 


43.16 


7 


(40,20,3) 








5.32e-3 


5.39e-3 


4.90 


12.8 


13.60 


5.83 


8 


(40,20,5) 








6.01e-3 


5.97e-3 


4.30 


12.1 


41.46 


17.92 


9 


(40,20,10) 








6.41e-3 


6.25e-3 


4.10 


8.90 


200.88 


60.21 



executing hardware floating-point arithmetic. For every example, we have gen- 
erated 10 random test cases as in the above. In executing the GPGCD method, 
we set u — 100 and a threshold of the 2-norm of the "update" vector in each 
iteration e = 1.0 x 10~ 8 ; in R_con_mulpoly, we set the tolerance e = 1.0 x 10~ 8 . 

Table [1] shows the results of the test. In each test, we have given several 
polynomials of the same degree as the input. The second column with (m, d, n) 
denotes the degree of input polynomials, degree of GCD, and the number of input 
polynomials, respectively. The columns with "STLN" are the data for the STLN- 
based method, while those with "GPGCD" are the data for the GPGCD method. 
"#Fail" is the number of "failed" cases such as: in the STLN-based method, the 
number of iterations exceeds 50 times (which is the built-in threshold in the 
program), while, in the GPGCD method, the number of iterations exceeds 100 
times. All the other data are the average over results for the "not failed" cases: 
"Error" is the sum of perturbation X)"=i ll^» — -^lll) where "ae — &" denotes 
a x 10~ b ; "^Iterations" is the number of iterations; "Time" is computing time 
in seconds. 

We see that, in the most of tests, both methods calculate approximate GCD 
with almost the same amount of perturbations. In the most of tests, the GPGCD 
method runs faster than STLN-based method. However, running time of the 
GPGCD method increases as much as that of the STLN-based method in some 
cases with relatively large number of iterations (such as Ex. 6). There is a case 
in which the GPGCD method does not converge (Ex. 6). Factors leading to such 
phenomena is under investigation. 

5 Concluding Remarks 

Based on our previous research ([17], [H]), we have extended our GPGCD 
method for more than two input polynomials with the real coefficients. We have 
shown that, at least theoretically, our algorithm properly calculates an approxi- 
mate GCD under certain conditions for multiple polynomial inputs. 
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Our experiments have shown that, in the case that the number of iterations is 
relatively small, the GPGCD method calculates an approximate GCD efficiently 
with almost the same amount of perturbations as the STLN-based method. 
However, computing time of the GPGCD method increases as the number of it- 
erations becomes larger; it suggests that we need to reduce the computing time 
of each iteration in the GPGCD method for the cases with relatively large num- 
ber of iterations. It is desirable to have more detailed experiments for analyzing 
stability, performance for input polynomials of larger degree, etc. 

For the future research, generalizing this result to polynomials with the com- 
plex coefficients will be among our next problems. It is also an interesting prob- 
lem how the choice of Pi affects the performance of the algorithm. Furthermore, 
one can also use arbitrary linear combination to transform gcd(Pi, P2, . . . , P n ) to 
gcd(Pi, a^Pi + • • • + a n P n ). This will reduce the size of the generalized Sylvester 
matrix and will be another approach for calculating approximate GCD. 
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